"""
Swift-Hohenberg PDE implementation

The Swift-Hohenberg equation is a model for pattern formation:
    ∂c/∂t = [ε - (kc² + ∇²)²]c + δc² - c³

This is a fourth-order PDE exhibiting rich pattern formation behavior.
"""

using ..Fields
using ..Grids

"""
    SwiftHohenbergPDE <: AbstractPDE

Swift-Hohenberg equation for pattern formation near onset.

The mathematical definition is:
    ∂c/∂t = [ε - (kc² + ∇²)²]c + δc² - c³

where c is the order parameter, ε is the bifurcation parameter (control parameter),
kc² is the critical wave number squared, and δ controls the quadratic nonlinearity.

# Parameters
- `rate`: Bifurcation parameter ε (default: 0.1)
- `kc2`: Critical wave number squared kc² (default: 1.0)  
- `delta`: Quadratic nonlinearity parameter δ (default: 1.0)
- `bc`: Boundary conditions for c (default: "auto_periodic_neumann")
- `bc_lap`: Boundary conditions for ∇²c (default: same as bc)

# Properties
- Fourth-order PDE with rich bifurcation structure
- ε > 0: Patterns emerge (supercritical bifurcation)
- ε < 0: Uniform state stable (subcritical regime)
- Exhibits hexagonal, stripe, and other patterns
- Energy functional: E = ∫[(ε - kc²)c² + c⁴/4 - δc³/3 + (∇²c)²/2] dV
"""
struct SwiftHohenbergPDE <: AbstractPDE
    rate::Float64        # ε - bifurcation parameter
    kc2::Float64         # kc² - critical wave number squared
    delta::Float64       # δ - quadratic nonlinearity
    bc::String           # boundary conditions for c
    bc_lap::String       # boundary conditions for ∇²c
    
    function SwiftHohenbergPDE(; rate=0.1, kc2=1.0, delta=1.0, 
                             bc="auto_periodic_neumann", bc_lap="")
        bc_lap_final = isempty(bc_lap) ? bc : bc_lap
        new(rate, kc2, delta, bc, bc_lap_final)
    end
end

"""
    evolution_rate(pde::SwiftHohenbergPDE, state::ScalarField, t::Float64)

Compute the evolution rate for the Swift-Hohenberg equation.

The equation is: ∂c/∂t = [ε - (kc² + ∇²)²]c + δc² - c³

Implementation:
1. Compute ∇²c
2. Compute (kc² + ∇²)²c = (kc²c + ∇²c)∇² = kc²∇²c + ∇⁴c  
3. Apply linear term: (ε - kc²)c - ∇⁴c - 2kc²∇²c
4. Add nonlinear terms: δc² - c³
"""
function evolution_rate(pde::SwiftHohenbergPDE, state::ScalarField, t::Float64)
    c = data(state)
    grid_obj = grid(state)
    
    # Step 1: Compute ∇²c
    c_laplace = laplacian(state)
    
    # Step 2: Compute ∇⁴c = ∇²(∇²c) (fourth-order term)
    c_biharmonic = laplacian(c_laplace)
    
    # Step 3: Linear part of the equation
    # (ε - kc²)c - 2kc²∇²c - ∇⁴c
    
    # Linear term: (ε - kc²)c
    linear_coeff = pde.rate - pde.kc2^2
    linear_term = ScalarField(linear_coeff * c, grid_obj)
    
    # Mixed term: -2kc²∇²c
    mixed_term = ScalarField(-2 * pde.kc2 * data(c_laplace), grid_obj)
    
    # Fourth-order term: -∇⁴c
    fourth_order_term = -c_biharmonic
    
    # Step 4: Nonlinear terms
    # Quadratic term: δc²
    quadratic_term = ScalarField(pde.delta * c.^2, grid_obj)
    
    # Cubic term: -c³
    cubic_term = ScalarField(-c.^3, grid_obj)
    
    # Step 5: Combine all terms
    result = linear_term + mixed_term + fourth_order_term + quadratic_term + cubic_term
    
    return result
end

"""
    free_energy(pde::SwiftHohenbergPDE, state::ScalarField)

Compute the free energy functional for the Swift-Hohenberg equation:
    E = ∫[f₀(c) + f₂(c, ∇²c)] dV
    
where:
    f₀(c) = (ε - kc²)c²/2 + δc³/3 + c⁴/4
    f₂(c, ∇²c) = (∇²c)²/2 + kc²|∇c|²
"""
function free_energy(pde::SwiftHohenbergPDE, state::ScalarField)
    c = data(state)
    
    # Local energy density f₀(c)
    local_energy = sum(
        (pde.rate - pde.kc2^2) * c.^2 / 2 + 
        pde.delta * c.^3 / 3 + 
        c.^4 / 4
    )
    
    # Gradient energy contributions
    c_laplace = laplacian(state)
    laplace_energy = sum(data(c_laplace).^2) / 2
    
    # |∇c|² term
    grad_c = gradient(state)
    grad_energy = 0.0
    if isa(grad_c, VectorField)
        grad_data = data(grad_c)
        if ndims(grad_data) >= 2
            for i in 1:size(grad_data, 1)
                grad_energy += sum(grad_data[i, :].^2)
            end
        else
            grad_energy += sum(grad_data.^2)
        end
    end
    grad_energy *= pde.kc2
    
    return local_energy + laplace_energy + grad_energy
end

"""
    is_pattern_forming(pde::SwiftHohenbergPDE)

Check if the parameters are in the pattern-forming regime (ε > 0).
"""
is_pattern_forming(pde::SwiftHohenbergPDE) = pde.rate > 0

"""
    get_critical_wavelength(pde::SwiftHohenbergPDE)

Get the critical wavelength λc = 2π/kc for pattern formation.
"""
get_critical_wavelength(pde::SwiftHohenbergPDE) = 2π / sqrt(pde.kc2)

"""
    get_linear_growth_rate(pde::SwiftHohenbergPDE, k)

Compute the linear growth rate for a mode with wave number k.
For the Swift-Hohenberg equation: σ(k) = ε - (k² - kc²)²
"""
function get_linear_growth_rate(pde::SwiftHohenbergPDE, k)
    return pde.rate - (k^2 - pde.kc2)^2
end

"""
    get_amplitude_equation_coefficients(pde::SwiftHohenbergPDE)

Get the coefficients for the amplitude equation near onset.
Returns (g₁, g₂) where the amplitude equation is: ∂A/∂T = εA + g₁|A|²A + g₂∇²A
"""
function get_amplitude_equation_coefficients(pde::SwiftHohenbergPDE)
    # These are derived from weakly nonlinear analysis
    g1 = -1.0  # Cubic coefficient (always stabilizing for SH)
    g2 = 1.0   # Spatial coupling in amplitude equation
    return (g1, g2)
end

"""
    is_stiff(pde::SwiftHohenbergPDE)

Swift-Hohenberg is stiff due to fourth-order derivatives.
"""
is_stiff(pde::SwiftHohenbergPDE) = true

"""
    requires_high_accuracy(pde::SwiftHohenbergPDE)

Pattern formation requires good accuracy to capture instabilities correctly.
"""
requires_high_accuracy(pde::SwiftHohenbergPDE) = abs(pde.rate) < 1.0

"""
    get_time_scale(pde::SwiftHohenbergPDE, grid_spacing)

Estimate the characteristic time scale for the Swift-Hohenberg equation.
"""
function get_time_scale(pde::SwiftHohenbergPDE, grid_spacing)
    # Time scale dominated by fourth-order term
    characteristic_time = grid_spacing^4
    
    # Modify based on proximity to bifurcation
    if abs(pde.rate) < 0.1
        characteristic_time /= abs(pde.rate)  # Slower near bifurcation
    end
    
    return characteristic_time
end

"""
    expression(pde::SwiftHohenbergPDE)

Return a string representation of the Swift-Hohenberg PDE.
"""
function expression(pde::SwiftHohenbergPDE)
    terms = String[]
    
    # Linear term: (ε - kc²)c
    linear_coeff = pde.rate - pde.kc2^2
    if linear_coeff != 0
        sign_str = linear_coeff >= 0 ? "" : "-"
        coeff_str = abs(linear_coeff) == 1.0 ? "" : "$(abs(linear_coeff))"
        push!(terms, "$(sign_str)$(coeff_str)c")
    end
    
    # Fourth-order and mixed terms: -∇²(2kc²c + ∇²c)
    if pde.kc2 != 0
        kc2_str = pde.kc2 == 1.0 ? "2" : "$(2*pde.kc2)"
        push!(terms, "-∇²($(kc2_str)c + ∇²c)")
    else
        push!(terms, "-∇⁴c")
    end
    
    # Quadratic term: δc²
    if pde.delta != 0
        delta_str = pde.delta == 1.0 ? "" : "$(pde.delta)"
        push!(terms, "$(delta_str)c²")
    end
    
    # Cubic term: -c³
    push!(terms, "-c³")
    
    return join(terms, " ")
end

# Export the PDE type and analysis functions
export SwiftHohenbergPDE, free_energy, is_pattern_forming, get_critical_wavelength,
       get_linear_growth_rate, get_amplitude_equation_coefficients, get_time_scale